      SUBROUTINE  JACROB (ESQ   ,HGT   ,TIME  ,USAT  ,USUN  ,RHO)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION       USAT(3)      ,USUN(3)
C
C.......................................................................
C
C  VERSION OF 20 JULY 1976
C
C   PURPOSE
C      JACROB FURNISHES ATMOSPHERIC DENSITY USING JACCHIA'S
C      1971 ATMOSPHERIC MODEL AS MODIFIED BY ROBERTS IN 1971.
C      COMPUTATIONS PERFORMED ARE TIME DEFINITION, EXOSPHERIC
C      TEMPERATURE, CORRECTIONS TO THE ATMOSPHERIC DENSITY, AND
C      ASSESSMENT AS TO WHICH HEIGHT RANGE TO USE IN COMPUTING THE
C      ATMOSPHERIC DENSITY.
C
C   REFERENCES
C      JACCHIA, L. J., REVISED STATIC MODELS OF THE THERMOSPHERE AND
C        EXOSPHERE WITH EMPIRICAL TEMPERATURE PROFILES, S.A.O., SPECIAL
C        REPORT NO. 332, MAY 5,1971.
C      ROBERTS, C. E., JR., AN ANALYTICAL MODEL FOR UPPER ATMOSPHERE
C        DENSITIES BASED UPON JACCHIA'S 1970 MODELS, CELESTIAL MECHANICS
C        4,(1971).
C
C    COMMON BLOCKS NEEDED
C       NONE
C
C    CALLING SEQUENCE EXPLANATION
C
C       VARIABLE       I/O     DESCRIPTION
C       --------       ---     -----------------------------------------
C       ESQ             I      ECCENTRICITY OF THE CENTRAL BODY SQUARED
C       HGT             I      SPACECRAFT HEIGHT (KM)
C       IERR            I      ERROR RETURN CODE FROM JACCWF
C       TIME            I      CURRENT TIME (FULL A.1 JULIAN DATE)
C       USAT(3)         I      UNIT POSITION VECTOR OF THE SPACECRAFT
C       USUN(3)         I      UNIT POSITION VECTOR OF THE SUN
C       RHO             O      ATMOSPHERIC DENSITY
C
C   JACROB IS CALLED FROM SUBROUTINE DRAG
C
C   SUBROUTINES AND FUNCTIONS REQUIRED
C     HIALT  - COMPUTES RHO AT ALTITUDES ABOVE 125 KM.
C     JACCWF - RETURNS GEOMAGNETIC INDEX AND EXOSPHERIC TEMPERATURE
C     LOWALT - COMPUTES RHO AT ALTITUDES BETWEEN 90 KM AND 125 KM
C
C.......................................................................
C
C                   INITIALIZE NECESSARY PARAMETERS
C
      DATA   DTR   /1.745329251994330D-2/
      DATA   TWOPI /6.283185307179586D0/
      DATA   PI    /3.141592653589793D0/
C
C                   MINIMUM HEIGHT IN KILOMETERS
      DATA   ZJ0   /90.D0/
C                   DENSITY AT MINIMUM HEIGHT IN GM/CM**3
      DATA   RHOZ  /3.46D-9/
C                   TEMPERATURE AT MINIMUM HEIGHT IN DEGREES KELVIN
      DATA   T0    /183.D0/
C
      IF (HGT .GE. ZJ0)  GO TO 10
C
C                   SPACECRAFT IS AT OR BELOW MINIMUM HEIGHT
      RHO = RHOZ
      GO TO 90
C
   10 IF (HGT .LT. 3000.D0)  GO TO 20
C
C                   THERE IS NO ATMOSPHERE ABOVE 3000 KM
C
      RHO = 0.D0
      GO TO 999
C
C                   CALCULATE SATELLITE GEODETIC LATITUDE
C
   20 D1  = DSQRT (USAT(1)**2 + USAT(2)**2)
      PHJ = DATAN (USAT(3)/(D1*(1.D0 - ESQ)))
C
C                   CALCULATE SUN HOUR ANGLE
C
      D2  = DSQRT (USUN(1)**2  +  USUN(2)**2)
      CX  = USUN(1) * USAT(1)  +  USUN(2) * USAT(2)
      TH  = DARCOS (CX/(D1*D2))
      HC  = USUN(1)*USAT(2)  -  USUN(2)*USAT(1)
      HCJ = TH / DTR  *  HC / DABS(HC)
C
C                  OBTAIN TEMPERATURE AND MAGNETIC VALUE AT CURRENT TIME
C     ...........
      CALL JACCWF (TIME   ,TEMP  ,MAG)
C     ...........
C
C                  COMPUTE THE DIURNAL EXOSPHERIC TEMPERATURE -- DL
C
      DEL0   = DATAN2 (USUN(3),D2)
      Q      = (HCJ + 43.D0) * DTR
      TAU    = (HCJ - 37.D0 + 6.D0   *DSIN(Q)) * DTR
      IF (TAU .LT. -PI)  TAU = TAU  +  TWOPI
      IF (TAU .GT.  PI)  TAU = TAU  -  TWOPI
      TAU    = .5D0 * TAU
      TH22   = ( DSIN (.5D0 * DABS (PHJ+DEL0) ) ) ** 2.2
      UN22   = ( DCOS (.5D0 * DABS (PHJ-DEL0) ) ) ** 2.2
      CSTAU3 = DCOS(TAU)**3
      DL     = TEMP * (1.D0 + .3D0 * (TH22 + (UN22 - TH22) * CSTAU3))
C
C                  GEOMAGNETIC  ACTIVITY
C
C   GEOMAGNETIC DISTURBANCES INFLUENCE THE ATMOSPHERIC DENSITY WITH A
C   TIME LAG OF 6.7 HOURS.  THE STORED VALUES ARE 10 TIMES THE ACTUAL
C   VALUES.
C
      TP1 = (MAG * 3  + 5) / 10
      TP1 = TP1 / 3.D0
C
C   THE GEOMAGNETIC TEMPERATURE CONTRIBUTION (DTG)  AND
C   THE BASE 10 LOGARITHM OF THE CONTRIBUTION TO RHO FROM GEOMAGNETIC
C   INFLUENCES ARE HANDLED DIFFERENTLY ABOVE AND BELOW 200KM.
C
      IF (HGT .LT. 200.D0)  GO TO 30
C
      DTG    = 28.D0 * TP1  +  .03D0 * DEXP(TP1)
      DL10RG = 0.D0
      GO TO 40
C
   30 DTG    = 14.D0  * TP1  +  .02D0  * DEXP(TP1)
      DL10RG = .012D0 * TP1  +  1.2D-5 * DEXP(TP1)
C
C   EXOSPHERIC TEMPERATURE IS THE TOTAL TEMPERATURE DUE TO DIURNAL AND
C   GEOMAGNETIC EFFECTS
C
   40 TINF = DL + DTG
C
C   SEMI-ANNUAL DENSITY VARIATION IS ADDED ON TO THE GEOMAGNETIC -- BOTH
C   IN BASE 10 LOGARITHMS
C
      PHI    = (TIME - 2436204.5D0) / 365.2422D0
      TSA    = TWOPI * PHI  +  6.035D0
      TSA    = PHI + .09544D0*((.5D0 + .5D0*DSIN(TSA))**1.65D0 - .5D0)
      B1     = TWOPI * TSA  +  4.137D0
      B2     = TWOPI * TSA * 2.D0  +  4.259D0
      GT     = .02835D0 + (.3817D0 + .17829D0*DSIN(B1)) * DSIN(B2)
      FZ     = (.5876D-6*HGT**2.331D0+.06328D0)*DEXP(-.002868D0*HGT)
      DL10RG = DL10RG  +  FZ * GT
C
C   SEASONAL LATITUDE VARIATION OF THE LOWER THERMOSPHERE IS ALSO ADDED
C   IN TO THE LOGARITHM OF DENSITY
C
      ZM90 = HGT - ZJ0
C
C                    HGT .GT. 440 GIVES UNDERFLOWS.
C
      IF (HGT .LT. 440.D0)  GO TO 50
C
      DL10RL = 0.D0
      GO TO 60
C
   50 EXPR = .014D0 * DEXP (-1.3D-3*ZM90**2) * DSIN (TWOPI*PHI+1.72D0)
      DL10RL = DSIN(PHJ)
      DL10RL = DL10RL * DABS(DL10RL) * EXPR * ZM90
C
   60 DL10RG = 10.D0 ** (DL10RG + DL10RL)
C
C                     INFLECTION POINT TEMPERATURE
C
      TX =371.6673D0+.0518806D0*TINF-294.3505D0*DEXP(TINF*(-.0021622D0))
C
C   THE OTHER CONTRIBUTIONS ARE CALCULATED BY SUBROUTINES HIALT OR
C   LOWALT, AS APPROPRIATE, AND ADDED INTO THE LOGARITHM OF RHO.
C
      IF (HGT .GT. 125.D0) GO TO 70
C
C          ......
      CALL LOWALT (HGT   ,RHOZ  ,TINF  ,TX    ,T0    ,ZJ0   ,RHO   )
C          ......
C
      GO TO 80
C
C          .....
   70 CALL HIALT (DEL0  ,HGT   ,PHJ   ,TINF  ,TX    ,T0    ,ZJ0   ,RHO)
C          .....
C
   80 RHO = RHO * DL10RG
C
C                      CONVERT RHO FROM CGS TO MKS UNITS.
C
   90 RHO = RHO*1.D12
C
C
  999 RETURN
      END
